function eipsp_find(varargin)
% eipsp_find is based on the template routine (below).
% finds the epsps or ipsps in a trace, and tags the times into the CONTROL database.
% 8/23/05 P. Manis

% basic template file for data analysis routines.
% gets the current datac window selection, and does the analysis
% for all selected record groups/blocks
% automatically refreshes the data arrays from the disk file if necessary
% and also shows how to generate a figure.
% 8/23/05 P. Manis

% the main routine is just a wrapper to handle the multiple selection and
% catch execution errors.
%
% Note: when testing, set the testflag to 1 so that we DO NOT do try/catch.
% you need to do this to catch the errors in the _go routine.
%
global CONTROL
sf = getmainselection; % read the list of selected traces

testflag = 1;

if(sf > 0)
    pflag = getplotflag; % read the plot flag (used in the routine template2)
    if(nargin > 0)
        cmd = varargin{1}; % read and pass the command
    else
        cmd = '';
    end;
    QueMessage('Template analysis routine', 1); % clear the que
    for i = 1:length(sf)
        if(testflag == 0)
            try
                eipsp_go(sf(i), pflag, cmd);
            catch
                fprintf(1, 'EIPSP Analysis Failed execution for sf = %d\n  ***[File: %s, records: %d - %d ]\n', ...
                    sf(i), CONTROL(sf(i)).filename, CONTROL(sf(i)).recbeg, CONTROL(sf(i)).recend);
            end;
        else
            eipsp_go(sf(i), 0, cmd); % also turn off plot flag when testing
        end;

    end;
end;



%--------------The real analysis routine------------------

function eipsp_go(sf, pflag, cmd)
% method of Vincent and Marty 1993 Neuron - compare two stretches of the
% trace and see when the difference crosses a threshold. When it does, mark
% that ("events"), and then continue until there is a peak (diff changes
% sign), mark that ("peaks"). Continue to the end of the trace marking both
% the onset of events and peaks. The parameters law and raw are the left
% average window and right average window, respectively, and may need to be
% adjusted depending on the rates of rise, the sampling rate, and the noise
% in the data. The gap is the space inbetween the two windows - can run
% from 0 to about 3 in different conditions. You have to play with the
% paramters to optimize the detection for a given set of
% recording conditions and traces.
% Aside from preallocating the result arrays (events, peaks), and computing
% the moving average by adding the new end point and subtracting the old
% first point in the average (which is much faster than computing the mean
% every time), there is little optimization of the code.
%

% access to all of the data and parameters
global DFILE ALLCH VOLTAGE CONTROL

[DFILE, DPAR, err] = ana_setup(DFILE, sf);
if(err ~= 0)
    return;
end;
cmdflag = 0;
switch(cmd)
    % set some control parameters
    case 'plot'
        nrecs = DPAR.records;
        cmdflag = 1;
    case 'try'
        nrecs = 1;
        cmdflag = 2;
    case 'reset'
        nrecs = 1;
        cmdflag = 3;
    otherwise
        nrecs = DPAR.records;
end;

% set the defaults
srcchan = [1,3];
law = 5;
raw = 5;
gap = 1;
thr = 0.1; % minimum event size (negative for IPSP detection)
peakslopefac = 0; % slope factor for peak detection
sthr = thr*peakslopefac;
t0 = 0;
tmax = 1000;
spike_thresh = 0;
lpf = 200; % low pass filter setting for EPSPs in current clamp, Hz.
notch = 1000/1.25; % try to notch some noise
minrt = 1; % min rise time in msec
minev = 0.2; % min peak size, in mv
if(cmdflag >= 2)
    CONTROL(sf).EIPSP=[];
end;
if(isfield(CONTROL, 'EIPSP') & ~isempty(CONTROL(sf).EIPSP))
    law = CONTROL(sf).EIPSP.leftwin;
    raw = CONTROL(sf).EIPSP.rightwin;
    gap = CONTROL(sf).EIPSP.gap;
    thr = CONTROL(sf).EIPSP.thr;
    peakslopefac = CONTROL(sf).EIPSP.peakslopefactor;
    if(isfield(CONTROL(sf).EIPSP, 'minev') & ~isempty(CONTROL(sf).EIPSP.minev)) % added later - may not be updated in structure yet
        minev = CONTROL(sf).EIPSP.minev;
    end;
    if(isfield(CONTROL(sf).EIPSP, 'minrt') & ~isempty(CONTROL(sf).EIPSP.minrt))
        minrt = CONTROL(sf).EIPSP.minrt;
    end;
    lpf = CONTROL(sf).EIPSP.lpf;
    if(isfield(CONTROL(sf).EIPSP, 'notch') & ~isempty(CONTROL(sf).EIPSP.notch))
        notch = CONTROL(sf).EIPSP.notch;
    end;
    srcchan = CONTROL(sf).EIPSP.srcchan;
end;
names = {'raw', 'law', 'gap', 'thr', 'peakslopefac', 'minev', 'minrt', ...
    'lpf', 'notch'};

for i = 1:length(names)
    eval([sprintf('if(length(%s) == 1), %s = %s*ones(1, length(srcchan)), end;', ...
        names{i}, names{i}, names{i})]);
end;
%for i = 1:length(names)
%    eval([sprintf('%s = reshape(%s, size(1, length(%s))', names{i})]);
%end;

% get user input if running single at a time.
if(length(getmainselection) == 1 & cmdflag ~= 1)
    prompt={'Source Channels', 'left window(pts)','right window (pts)', 'gap (pts)', 'threshold (slope)', ...
        'peakslopefac', 'min event size (mV)', 'min rise time (msec)', ...
        'LPF (Hz)', 'Notch freq (Hz)'};
    def={num2str(srcchan), num2str(law), num2str(raw), num2str(gap), num2str(thr), ...
        num2str(peakslopefac), num2str(minev), num2str(minrt), ...
        num2str(lpf), num2str(notch)};
    dlgTitle='Parameters for EIPSP-find function';
    lineNo=1;
    ans=inputdlg(prompt,dlgTitle,lineNo,def);
    if(~isempty(ans))
        srcchan = str2num(ans{1});
        law = str2num(ans{2});
        raw = str2num(ans{3});
        gap = str2num(ans{4});
        thr = str2num(ans{5});
        peakslopefac = str2num(ans{6});
        minev = str2num(ans{7});
        minrt = str2num(ans{8});
        lpf = str2num(ans{9});
        notch = str2num(ans{10});
    else
        return;
    end;
end;



if(cmdflag ~= 1) % default for full analysis

    ev=cell(length(srcchan), DPAR.records);
    tev = ev;
    pk = ev;
    tpk = ev;

    jsrc = 0;
    for isrc = srcchan; % for all of the channels that need to be analyzed this way...
    jsrc = jsrc + 1;
        % this next routine reads all of the data and information we need,
        % including calculating some standard stuff that is put into DPAR.
        % you will want DPAR.time especially - it is the time base that goes
        % with the data.

        btr{isrc}  = ALLCH{isrc};
        trace = ALLCH{isrc};
        for i=1 : nrecs
            fsamp = 1000/DPAR.RATES(i); % get sampling frequency
            fco1 = 1000;		% highpass cutoff frequency in Hz
            fco2 = lpf(jsrc); % low-pass cutoff frequency in Hz
            wco1 = fco1/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
            wco2 = fco2/(fsamp/2);
            fn1 = notch(jsrc)+0.02*notch(jsrc);
            fn2 = notch(jsrc)-0.02*notch(jsrc);
            nco1 = fn1/(fsamp/2);
            nco2 = fn2/(fsamp/2);
            if(wco1 < 1 & wco2 < 1 & nco1 < 1 & nco2 < 1) % if wco is > 1 then this is not a filter!
                [bn, an] = fir1(24, [nco2 nco1], 'stop');
                tracel = filtfilt(bn, an, trace(i,:)-mean(trace(i,:)));
                [b, a] = fir1(24, wco2); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
                tracel = filtfilt(b, a, tracel); % filter all the traces...
                trace(i,:) = tracel;
                btr{isrc}(i,:) = tracel; % save smoothed trace...
            else
                fprintf(1, 'ng: wco1 = %f wco2 = %f\n', wco1, wco2);
            end


            events = NaN*zeros(1000, 1);
            peaks = events;
            te = events;
            tp = events;
            nevent = 0;
            maxj = size(trace,2)-(law+raw+gap);
            la = 0;
            ra = 0;
            first = 0;
            thistr = trace(i,:); % make local allocation of the trace
            first = 0;
            fprintf(1, 'Ch: %d  Rec: %d\n', isrc, i);
            j = law(jsrc)+1;
            while (j <= maxj)
                jrangel = [j-law(jsrc):j];
                jranger = [j+gap(jsrc):j+gap(jsrc)+raw(jsrc)];
                j = j + 1;
                if(first == 0)
                    la = mean(thistr(jrangel));
                    ra = mean(thistr(jranger));
                    first = 1;
                else
                    la = la + (thistr(max(jrangel)+1) - thistr(min(jrangel)))/law(jsrc);
                    jrangel = [jrangel(2:end) max(jrangel)+1];
                    ra = ra + (thistr(max(jranger)+1) - thistr(min(jranger)))/raw(jsrc);
                    jranger = [jranger(2:end) max(jranger)+1];
                end;
                if(((la - ra) > thr(jsrc)))
                    nevent = nevent + 1;
                    events(nevent) = thistr(j); % mark start of event
                    peaks(nevent) = NaN; % not found yet
                    te(nevent) = DPAR.time(i,j);
                    tp(nevent) = NaN;
                    while((((la - ra) > sthr) & j < maxj)) % wait for it to switch sign.
                        la = la + (thistr(max(jrangel)+1) - thistr(min(jrangel)))/law(jsrc);
                        ra = ra + (thistr(max(jranger)+1) - thistr(min(jranger)))/raw(jsrc);
                        j = j + 1; % advance through the trace
                        jrangel = [j-law(jsrc):j];
                        jranger = [j+gap(jsrc):j+gap(jsrc)+raw(jsrc)];
                        first = 0;
                    end;
                    peaks(nevent) = thistr(j);
                    tp(nevent) = DPAR.time(i,j);
                    if((abs(peaks(nevent)-events(nevent)) < minev(jsrc)) | ...
                            (tp(nevent)-te(nevent)) < minrt(jsrc)) % filter small or too-short events
                        peaks(nevent) = NaN;
                        events(nevent) = NaN;
                        te(nevent) = NaN;
                        tp(nevent) = NaN;
                        nevent = nevent - 1;
                    end;
                end;

            end;
            [iev jev] = find(~isnan(events)); % reduce to just the captured events  & real stuff
            ev{isrc,i} = events(iev);
            tev{isrc,i} = te(iev);
            pk{isrc,i} = peaks(iev); % may include NaN's if peak fell after the end of the trace...
            tpk{isrc,i} = tp(iev);
        end;
    end;
else % end of empty command... get the data from the CONTROL database and then plot it.
    jsrc = 0;
    for isrc = srcchan
        jsrc = jsrc + 1;
        a = eval([sprintf('CONTROL(sf).EIPSP.ev%d;', isrc)]);
        ta = eval([sprintf('CONTROL(sf).EIPSP.tev%d;', isrc)]);
        b = eval([sprintf('CONTROL(sf).EIPSP.pk%d;', isrc)]);
        tb = eval([sprintf('CONTROL(sf).EIPSP.tpk%d;', isrc)]);
        for k = 1:length(a)
            ev{isrc, k} = [a{k}];
            tev{isrc, k} = [ta{k}];
            pk{isrc, k} = [b{k}];
            tpk{isrc,k} = [tb{k}];
        end;
        % filter the data
        btr{isrc}  = ALLCH{isrc};
        trace = ALLCH{isrc};
        for i=1 : nrecs
            fsamp = 1000/DPAR.RATES(i); % get sampling frequency
            fco1 = 1000;		% high-pass cutoff frequency in Hz
            fco2 = lpf(jsrc); % low-pass cutoff frequency in Hz
            wco1 = fco1/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
            wco2 = fco2/(fsamp/2);

            if(wco1 < 1 & wco2 < 1) % if wco is > 1 then this is not a filter!
                [b, a] = fir1(8, wco2); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
                tracel = filtfilt(b, a, trace(i,:)-mean(trace(i,:))); % filter all the traces...
                trace(i,:) = tracel;
                btr{isrc}(i,:) = tracel; % save smoothed trace...
            else
                fprintf(1, 'ng: wco1 = %f wco2 = %f\n', wco1, wco2);
            end
        end;
    end;
end;

% set up a figure window (or you could grab an existing one)
hf = findobj('tag', 'EIPSP'); % find the tagged window with it's name
if(isempty(hf)) % doesn't exist
    hf = figure();
    set(hf, 'tag', 'EIPSP');
    set(hf, 'NumberTitle', 'off');
    set(hf, 'Name', 'EIPSP');
end;
figure(hf);


clf;% clear it.


% Define a text area and print some standard junk.
subplot('position', [0.1, 0.90, 0.50, 0.095]);
axis([0,1,0,1])
axis('off')
ht(1)=text(0,0.80,sprintf('%-12s R[%d:%d]     %-8s  [%s]',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol, date), 'Fontsize', 10);
set(ht(1), 'interpreter', 'none'); % un-TeX the line - this is a filename and won't have tex chars, but might have an underscore.
text(0,0.6,sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 8);
text(0,0.4,sprintf('Ihold:%6.2f %s    RMP: %6.2f %s, Rin: %8.3f M\\Omega', ...
    CONTROL(sf).iHold,CONTROL(sf).I_Unit, CONTROL(sf).Rmp, CONTROL(sf).V_Unit, CONTROL(sf).Rin), 'FontSize', 8);


% define some subplots and plot some data

subplot('position', [0.1, 0.48, 0.50, 0.320]);
tisih = [0:0.5:50];
isi1 = [];
for i = 1:nrecs
    plot(DPAR.time(i,:), btr{srcchan(1)}(i,:)); % raw, unfiltered trace...?
    hold on
    plot(tev{srcchan(1),i}, ev{srcchan(1),i}, 'rx');
    plot(tpk{srcchan(1),i}, pk{srcchan(1),i}, 'bs');
    isi1 = [isi1 diff([tev{srcchan(1), i}]')];
end;
h1 = gca; % get a handle to the axes.
% label the plot and stuff if you need
xlabel('Time (ms)', 'fontsize', 9);
ylabel(sprintf('CH%d V (mV)', srcchan(1)), 'fontsize', 9);
set(gca, 'fontsize', 8);
box off % turn off the outer right and top axes marks.

subplot('position', [0.7, 0.72, 0.27, 0.22]);
hisi = histc(isi1, tisih);
if(~isempty(isi1))
    bar(tisih, hisi);
    set(gca, 'Xlim', [0 50]);
    ylabel(sprintf('CH %d ISI', srcchan(1)),'fontsize', 9);
    set(gca, 'fontsize', 8);
else
    text(0.5, 0.5, 'No events', 'horizontalalignment', 'center');
end;

if(length(srcchan) > 1)
    subplot('position', [0.1, 0.0750, 0.50, 0.320]);
    isi2 = [];
    for i = 1:nrecs
        plot(DPAR.time(i,:), btr{srcchan(2)}(i,:)); % raw, unfiltered trace...?
        hold on
        plot(tev{srcchan(2),i}, ev{srcchan(2),i}, 'rx');
        plot(tpk{srcchan(2),i}, pk{srcchan(2),i}, 'bs');
        isi2 = [isi2 diff([tev{srcchan(2), i}]')];
    end;

    h1 = gca; % get a handle to the axes.
    % label the plot and stuff if you need
    xlabel('Time (ms)', 'fontsize', 9);
    ylabel(sprintf('CH%d V (mV)', srcchan(2)), 'fontsize', 9);
    box off % turn off the outer right and top axes marks.
    set(gca, 'fontsize', 8);

    subplot('position', [0.7, 0.480, 0.27, 0.22]);
    if(~isempty(isi2))
        hisi2 = histc(isi2, tisih);
        bar(tisih, hisi2);
        set(gca, 'Xlim', [0 50]);
        ylabel(sprintf('CH %d ISI', srcchan(2)), 'fontsize', 9);
        set(gca, 'fontsize', 8);
    else
        text(0.5, 0.5, 'No events', 'horizontalalignment', 'center');
    end;
end;

if(cmdflag ~= 1)
    CONTROL(sf).EIPSP=[];
    CONTROL(sf).EIPSP.leftwin = law;
    CONTROL(sf).EIPSP.rightwin = raw;
    CONTROL(sf).EIPSP.gap = gap;
    CONTROL(sf).EIPSP.thr = thr;
    CONTROL(sf).EIPSP.peakslopefactor = peakslopefac;
    CONTROL(sf).EIPSP.minev = minev;
    CONTROL(sf).EIPSP.lpf = lpf;
    CONTROL(sf).EIPSP.srcchan = srcchan;
    for isrc = srcchan
        eval([sprintf('CONTROL(sf).EIPSP.tev%d = {tev{isrc,:}}',isrc)]);
        eval([sprintf('CONTROL(sf).EIPSP.ev%d = {ev{isrc,:}}',isrc)]);
        eval([sprintf('CONTROL(sf).EIPSP.tpk%d = {tpk{isrc,:}}',isrc)]);
        eval([sprintf('CONTROL(sf).EIPSP.pk%d = {pk{isrc,:}}',isrc)]);
    end;
end;

if(length(srcchan) == 2)
    t1 = cell(nrecs,1);
    t3 = cell(nrecs,1);
    for i = 1:nrecs
        t = [tev{1,i}];
        t1{i} = t;
        t  = [tev{3,i}];
        t3{i} = t;
    end;
    h3 = subplot('Position', [0.7, 0.075, 0.27, 0.320]);
    simxcorr2(t1, t3, h3);
    ylabel(sprintf('XCorr ch: %d-%d', srcchan(1), srcchan(2)), 'fontsize', 9);
    xlabel('T (ms)', 'fontsize', 9);
set(gca, 'fontsize', 8);
end;
if(pflag) % plotflag and enable figure printing
    orient landscape % make the graph landscape on the paper.
    print;
end;
